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SCIENCE 


We  present  a  stable  method  for  the  inverse  scattering  problem  of  the  Helmholtz 
equation  in  two  dimensions.  The  algorithm  requires  single-frequency  scattering 
data,  and  is  an  iterative  procedure  which  resembles  the  process  of  layer- stripping. 
The  inversion  method  is  based  on  the  observation  that  the  iU-posedness  of  the 
inverse  scattering  problem  causes  it  to  be  almost  linear  in  certain  regimes.  In 
these  regimes,  the  algorithm  solves  the  resulting  quasi-linear  equations  to  produce 
approximate  solution  to  the  inverse  problem  within  a  narrow  circular  layer  sur¬ 
rounding  the  yet  unrecovered  part  of  the  scatterer.  This  approximation  is  used  to 
linearize  the  underlying  narrow  circular  strip;  in  the  process,  the  previously  ob¬ 
tained  solution  is  refined.  The  performance  of  the  algorithm  is  demonstrated  with 
several  numerical  examples  for  the  special  case  of  radially  symmetric  scatterers. 
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1  Introduction 


In  [7],  a  stable  method  is  presented  that  solves  the  inverse  scattering  problem  for 
the  Helmholtz  equation  in  two  dimensions.  The  algorithm  is  an  iterative  proce¬ 
dure,  and  requires  multi-frequency  scattering  data.  In  this  paper,  we  propose  and 
test  a  single-frequency  inversion  method  that  solves  the  same  inverse  scattering 
problem. 

There  are  two  major  difficulties  in  the  solution  of  an  inverse  problem  naturally 
formulated  as  a  system  of  nonlinear  equations:  ill-posedness  and  local  minima. 
They  seem  to  act  indissolubly  together  to  make  the  investigation  of  the  inverse 
problem  even  more  difficult.  Here  is  a  typical  situation  where  we  are  confronted 
with  a  dilemma.  On  one  hand,  the  lower  the  frequency  employed  in  the  inversion, 
the  fewer  the  local  minima  there  are  in  the  nonlinear  system,  but  the  more  ill- 
posed  the  problem  becomes.  On  the  other  hand,  to  reduce  the  ill-posedness  and 
therefore  to  increase  the  resolution  of  the  inversion,  a  relatively  high  frequency  is 
desired;  but  the  higher  the  frequency,  the  more  numerous  the  local  minima  there 
are  in  the  nonlinear  system. 

When  high  resolution  is  sought  in  a  reconstruction,  it  can  only  be  achieved 
with  high-frequency.  Therefore,  the  problem  of  local  minima  must  be  dealt  with 
there.  The  fact  that  no  optimization-based  techniques  have  so  far  been  made 
reliable  shows  that  it  is  difficult  to  directly  attack  the  problem.  One  of  the 
well-known  and  systematic  elforts  made  in  solving  the  problem  of  local  minima 
is  the  method  of  layer-stripping  (see,  for  example,  see  [1],  [2],  [3]),  in  both  the 
time  and  frequency  domains.  In  this  approach,  the  inverse,  nonlinear  problem 
is  often  reformulated  as  a  nonlinear  ordinary  differential  equation,  usually  of 
Riccati  type  in  the  frequency  domain,  for  the  scattering  matrix  in  the  context 
of  inverse  scattering  (see  [6]),  or  for  the  Dirichlet-to-Neumann  mapping  in  the 
context  of  electrical  impedance  imaging  (see  [1],  [3]).  The  measured  data  become 
the  initial  values  of  the  ODE.  The  initial  value  problem  is  solved,  which  can 
be  interpreted  as  layer-stripping  the  object  being  imaged,  and  the  problem  of 
local  minima  is  simply  eliminated.  However,  the  propagation  of  the  measured 
information  into  the  object  is  unstable,  for  the  initial  value  problem  solved  here 
is  similar  and  connected  to  the  Cauchy  problem  of  the  Laplace  equation  in  the 
case  of  electrical  impedance  imaging,  and  of  the  Helmholtz  equation  in  the  case 
of  inverse  scattering,  both  being  ill-posed. 

The  single-frequency  inversion  method  to  be  presented  in  this  paper,  as  well 
as  the  multi-frequency  method  of  [7],  is  based  on  the  observation  that  the  ill- 
posedness  of  the  inverse  scattering  problem  causes  it  to  be  almost  linear  in  cer¬ 
tain  regimes.  In  these  regimes,  the  algorithm  solves  the  resulting  quasi-linear 
equations  to  produce  approximate  solution  to  the  inverse  problem  within  a  nar¬ 
row  circular  layer  surrounding  the  yet  unrecovered  part  of  the  scatterer.  This 
approximation  is  used  to  linearize  the  underlying  narrow  circular  strip;  in  the 
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process,  the  previously  obtained  solution  is  refined. 

The  underlying  principle  of  physics  used  here  is  referred  to  as  the  skin  effect. 
For  elliptic  partial  differential  equations  in  two  (or  higher)  dimensions,  a  solution 
highly  oscillating  in  a  spatial  direction  decays  or  grows  rapidly  (usually  expo¬ 
nentially)  in  directions  perpendicular  to  it.  This  is  the  case  for  any  oscillatory 
solution  of  the  Laplace  equation,  or  for  an  evanescent  wave  of  the  Helmholtz 
equation.  Therefore,  when  the  object  is  probed  on  its  boundary  with  a  wave 
highly  oscillating  across  its  circumference,  only  a  thin  layer  of  the  object  is  pene¬ 
trated,  hence  the  term  skin  effect.  Corresponding  to  this  decaying  incident  field, 
the  scattered  field  which  we  measure  on  the  boundary  contains  information  of  the 
object  in  that  thin  layer.  Such  a  measurement  is  entirely  inadequate  to  determine 
the  whole  object — this  obviously  is  a  severely  ill-posed  problem.  But  the  mea¬ 
surement  may  be  used  to  approximately  determine  the  object  just  in  that  thin 
layer.  We  can  in  fact  do  so  since  the  problem  of  obtaining  the  thin  layer  from  the 
measurement  is  essentially  linear.  That  is  how  the  recursive  linearization  process 
starts. 

Less  highly  oscillatory  incident  waves  are  then  irradiated  onto  the  object. 
While  the  probing  energy  penetrates  a  thicker  layer  of  the  object,  the  relationship 
between  the  measurements  and  the  parameters  to  be  recovered  in  the  deeper 
layer  becomes  more  nonlinear.  These  nonlinear  equations  can  be  considered  as 
perturbations  to  the  already  solved  equations  at  the  previous  layers,  and  therefore 
can  be  continually  and  recursively  linearized  with  the  standard  techniques  of 
regular  perturbations.  Thus,  the  recursive  linearization  is  a  continuation  method 
on  a  parameter  of  the  the  incident  field  (see  (8)  and  the  integer  m  there)  which 
controls  the  depth  of  its  penetration. 

The  inversion  method  is  inherently  multi-dimensional  in  the  sense  that  in  one¬ 
dimensional  space  the  Cauchy  problem  for  an  elliptic  equation  is  not  ill-posed, 
that  there  is  no  decaying  mode  or  skin  effect  there,  and  therefore  that  the  method 
has  no  one-dimensional  version.  Also,  had  the  problem  not  been  ill-posed  in  two 
and  higher  dimensions,  it  might  have  been  more  difficult  to  design  an  algorithm 
to  solve  the  inverse  problem  with  numerous  local  minima. 

The  main  purpose  of  this  paper  is  to  propose  a  single-frequency  inversion 
method,  and  to  demonstrate  its  preliminary  numerical  results.  The  plan  of  the 
paper  is  as  follows:  in  Section  2,  we  reformulate  the  inverse  scattering  problem 
using  the  concept  of  illuminated  areas;  in  Section  3,  we  describe  the  inversion 
method  for  the  Helmholtz  equation  in  two  dimensions.  In  Section  4,  we  im¬ 
plement  the  method  numerically  for  the  special  case  of  cylindrically  symmetric 
scatterers,  and  present  the  numerical  results.  Finally,  in  Section  5,  we  conclude 
with  brief  remarks  on  the  performance  of  the  inversion  algorithm  and  discuss  its 
generalizations. 

Remark  1.1  Although  our  numerical  experiments  demonstrate  convergence  and 
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stability  of  the  inversion  algorithm,  its  analysis  is  presently  incomplete.  There¬ 
fore,  the  results  presented  in  this  paper  should  be  viewed  as  experimental. 

2  Skin  Effect  and  Illuminated  Areas 

In  this  section,  we  formulate  the  forward  and  inverse  scattering  problems  for  the 
Helmholtz  equation 

A(l>{x)  +  k^{l q{x))<f>{x)  =  0  (1) 

in  two  dimensions,  and  describe  the  skin  effect  and  related  concepts. 


2.1  The  Scattering  Problem 

In  (1),  we  assume  that  Ar  is  a  positive  number,  and  ^  is  a  smooth  function  with 
compact  support  f2  C  we  will  be  referring  to  the  function  q  as  the  scatterer, 
or  the  forward  model.  We  will  be  considering  solutions  of  (1)  of  the  form 

<l>{x)  =  <f>o{x)  +  ik{x),  (2) 

with  <j>Q  the  incident,  or  the  incoming,  field  satisfying  in  f)  the  homogeneous 
Helmholtz  equation 

A<j)o{x)  +  k^4>o{x)  =  0,  (3) 

and  with  tp  the  scattered  field  subject  to  the  outgoing  (Sommerfeld)  radiation 
condition 

We  will  be  referring  to  the  determination  of  the  scattered  field  from  a  given 
incoming  field  as  a  forward  scattering  problem.  It  is  well-known  (see,  for  example, 
[4])  that  the  forward  scattering  problem  is  well-posed. 

It  is  also  well-known  that  the  scattering  problem  can  be  reformulated  as  the 
Lippmann- Schwinger  equation  for  the  scattered  field  ip 

Hx)  =  -k^  I  (MOiMO  +  mn.  (5) 

J  Q 

where  Gk  is  the  free-space  Green’s  function  for  the  Helmholtz  equation.  The 
inverse  scattering  problem  is  to  determine  the  scatterer  q  inside  the  domain  Cl 
from  measurements  of  the  scattered  fields  outside  the  scatterer.  For  definiteness, 
we  assume  that  the  measurements  are  taken  on  the  boundary  of  Cl: 

{  'Pix)\x^dci  }.  (6) 

The  inverse  problem  is  nonlinear  since  both  the  scatterer  and  the  scattered  field 
is  unknown  inside  fi,  and  there  is  a  product  of  them  in  the  integral  equation  (5). 
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2.2  Skin  Effect  and  Weak  Scattering 

For  simplicity  and  without  loss  of  generality,  we  assume  that  fl  C  is  a  disk 
D{w)  defined  by  the  formula 


=  {  (x,J/),  ||(x,t/)l|  <  G7  },  (7) 


for  some  ru  >  0.  As  is  well-known,  inside  D{w),  the  solutions  of  (3)  are  of  the 
form 


=  Jm{kr)  ■ 


(8) 


where  m  is  integer,  referred  to  as  the  propagation  number  hereafter;  is  the 
Bessel  function  of  order  m,  and  =  (— 1)™  Jm-  Furthermore,  any  incident  field 

to  D{uj)  is  a  linear  combination  of  (8);  therefore,  we  only  need  to  consider  (8)  as 
the  incident  fields  for  the  purpose  of  inverse  scattering.  For  n  greater  than  ifccj, 
the  Bessel  function  Jn{kr)  decays  monotonically  as  the  radius  r  decreases  from 
ZD  to  zero.  In  fact. 


Jn{kr) 


1 

. . . — 

^/2'irn 


(9) 


for  kr  substantially  smaller  than  n.  Therefore,  for  |m|  substantially  greater  than 
kw,  the  incident  field  (8)  decays  rapidly  in  the  neighborhood  of  r  =  tu.  It  can 
only  penetrate  the  skin  of  the  scatterer  which  interacts  with  the  incident  field, 
and  produces  a  weak  scattered  field  Thus,  the  Born  approximation  can  be 
used  to  linearize  the  relationship  between  q  and  xj^.  In  other  words,  dropping  the 
term  containing  the  product  of  q  and  ip  in  the  Lippmann-Schwinger  equation  (5), 
to  the  second  order  of  ip^'^\  we  have 


=  -e  j  (10) 

for  \m\  substantially  greater  than  kzj.  Note  that  the  major  contribution  of  the 
integral  over  D{zd)  comes  from  a  thin  layer  of  the  disk  near  the  boundary. 


2.3  Illuminated  Areas 

As  is  shown  in  the  preceding  section,  the  incident  field  illuminates  a  thin 
layer  of  the  scatterer  for  |m|  substantially  greater  than  kw.  When  \m\  decreases, 
the  layer  the  incident  field  penetrates  becomes  thicker.  In  this  section,  we  describe 
this  process  and  introduce  necessary  notation. 

As  is  well-known  (see,  for  example,  [5]),  for  an  integer  n  >  0,  Jn{z)  is  an 
oscillatory  function  of  2  for  z  >  n.  It  decays  as  z  decreases  in  [0,n],  and  it 
decays  to  zero  at  the  rate  z"  as  z  goes  to  zero.  Therefore,  given  a  small  number 
e  >  0,  there  exists  z„  >  0  such  that  \Jn{z)\  <  t  for  all  z  <  Zy^.  Furthermore,  z„ 
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is  roughly  of  the  order  n,  and  decreases  as  n  decreases.  Therefore,  in  the  disk 
D{w),  the  incident  field 

4"’(i',e)  =  (11) 

illuminates  an  annular  area 

Ao(n,  A:)  =  {  (r,  I  r„  <  r  <  07  },  (12) 

where  k  ■  rn  =  Zn^  and  r„  decreases  as  n  decreases.  Figure  1  shows  sixteen 


Figure  1:  Absolute  Values  of  the  Incident  Fields. 

curves  depicting  n  +  |<^cr^l  for  n  =  0, 1, . . . ,  15,  fc  =  6.2,  as  functions  of  the  radius 
r  €  [0,  tt].  For  n  >  1,  the  flat  portion  of  each  curve  corresponds  to  the  area  that 
is  near  the  center  of  the  scatterer  and  is  not  illuminated  by  the  incident  field 
The  area  that  is  not  illuminated  by  is  essentially  the  center  r  =  0,  whereas 
illuminates  the  entire  scatterer. 

Remark  2.1  It  is  possible  to  rigorously  define  and  calculate  r„,  a  process  omitted 
here  since  it  is  tedious  and  its  result  is  irrelevant  to  the  subsequent  discussions. 

Definition  2.2  Let  V  be  a  domain  in  R.^,  and  g  :  V  ^  ^  a  smooth  function. 
For  e  >  0,  we  refer  to 

{  (a:,y),  |^(x,j/)|  >  e  }  C  V  (13) 
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as  the  e-essential  support  of  g.  We  will  omit  e  from  the  expression  “e-essential” 
if  e  is  a  small  number  and  its  explicit  reference  there  is  superfluous. 


It  follows  from  Definition  2.2  that  the  area  Ao{n,  k)  illuminated  by  the  incident 
field  <l>^^  is  the  essential  support  of  the  incident  field  in  D{w).  We  similarly 
define  the  illuminated  eirea  of  the  total  field 

+  (14) 

as  the  essential  support  of  the  total  field  in  D(zu),  and  denote  it  by  A{m,k). 


Figure  2:  Absolute  Values  of  the  total  Fields. 


Therefore, 

A{n,k)=  IJ  A{m,k)  (15) 

\m\>n 

is  the  illuminated  axea  of  the  scatterer  by  the  total  fields 

{  |m|  >  n  }.  (16) 

We  also  refer  to  A(n,  k)  as  the  illuminated  area  by  the  scattering  experiment 

with  incident  fields  (16).  Obviously,  D(cu)  D  A(0,  k)  D  A(l,  k)  D  A{2,  k)  D 

Figure  2  shows  sixteen  curves  depicting  n  +  for  n  =  0,1,...,  15  as 
functions  of  the  radius  r,  where  is  the  total  field  corresponding  to  the  forward 
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scattering  problem  with  w  =  tt,  k  =  Q.2  and  a  cylindrically  symmetric  scatterer 
q  defined  by  the  formula 

q(^r)  =  0.26  •  [0.55cos(4r)  —  0.44sin(8r)  +  0.25sin(12r)  +  0.14cos(16r)].  (17) 

For  n  >  1,  the  flat  portion  of  each  curve  corresponds  to  the  area  that  is  near 
the  center  of  the  scatterer  and  is  not  illuminated  by  the  total  field  The 

area  that  is  not  illuminated  by  is  essentially  the  center  r  =  0,  whereas 
illuminates  the  entire  scatterer. 

2.4  Reformulating  the  Inverse  Scattering  Problem 

The  essence  of  inversion  is  to  determine  a  forward  model  q,  which  may  not  be 
the  same  as  q,  but  which  produces,  exactly  or  approximately,  the  scattering 
measurement 

{  ‘tp^‘^\x)\a;edD(,ka)  }  (18) 

in  the  scattering  experiment  with  the  incident  field  <l>^^  for  all  integer  m.  More¬ 
over,  for  the  sake  of  stability,  we  require  that  q  has  small  norm.  In  this  section, 
we  formalize  these  two  requirements. 

Definition  2.3  Suppose  that  in  a  scattering  experiment,  the  two  forward  models 
g  :  I— K,  and  9  :  fl  »->  R,  produce  respectively  the  two  scattering  measurements 
ipldQ  and  For  a  given  e  >  0,  the  two  scattering  measurements  are  said  to 

be  e-essentially  the  same  if 

IKV’lafi  - '*^|an)||2  <  £•  (19) 

We  will  omit  e  from  the  expression  “e-essentially”  if  e  is  meant  to  be  a  small 
number  and  its  explicit  reference  there  is  superfluous. 

Definition  2.4  To  a  scattering  experiment,  the  two  forward  models  q  and  q  are 
said  to  look  (essentially)  the  same  if  they  produce  essentially  the  same  scattering 
measurements  in  the  experiment. 

Definition  2.5  A  forward  model  q  is  said  to  be  observable,  or  an  observable  part 
of  the  original  scatterer  q,  to  a  scattering  experiment  at  frequency  k,  if  it  looks 
the  same  as  the  original  q,  and  its  norm  is  the  least  among  those  that  look  the 
same  as  q. 

It  follows  from  the  Lippmann-Schwinger  equation  (5)  and  Definition  2.5  that  an 
inversion  using  the  incident  fields 

{  4’">,  |m|  >  n  }  (20) 

can  only  at  best  produce  observable  part  of  q  in  A{n,  k). 
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Definition  2.6  Given  an  integer  n  >0,  we  denote  by  E{n)  the  scattering  exper¬ 
iment  with  the  incident  fields  (20),  and  by  qn  :  A{n,k)  the  observable  part 

of  q  corresponding  to  the  experiment  E{n)  (here  it  is  assumed  that  the  observable 
forward  model  qn  exists  and  is  unique). 

With  qn  and  A{n,  k),  the  scattering  problem  (5)  can  be  modified  as  (see  Remark 

2.10) 

f  G,(x,  (21) 

J  A(n^k) 

for  |m|  >  n.  And  therefore  the  inverse  problem  is  to  obtain  which  satisfies 

/  G,(x,  jli? ,  (22) 

^A(0,k) 

with  prescribed  scattering  data 

{  U€9r>(tz7)?  ^  integer  }.  (23) 

Remark  2.7  More  precisely ^  go  has  the  least  norm  among  those  forward  mod¬ 
els  g  whose  scattering  data 

{  }  (24) 

satisfy 

(E  li(l(’'”’l8C(-)  -  1*‘’"’|90(-))II?)  '  <  (25) 

for  a  prescribed  small  number  e  >  0. 

Remark  2.8  Suppose  that  v  \  is  a  positive  number,  and  that  the  forward 
model  q  is  observable  to  the  scattering  experiment  with  the  incoming  fields 

{  ^01,  4>o2  },  (26) 

whereas  the  forward  model  q  is  observable  to  the  scattering  experiment  with  the 
incoming  fields 

{  <f>Qli  •  <f>Q2  }-  (27) 

Then,  the  two  forward  models  are  not  the  same  in  general,  because  the  cost  func¬ 
tions  (25)  corresponding  to  the  two  sets  of  incident  fields  are  different. 

Remark  2.9  In  general,  qn  :  A(n,  k)  R,^  and  :  i4(n  —  1,  A:)  are  not 

identical  in  A{n  —  1,  k),  but  both  look  the  same  to  the  scattering  experiment  E{n). 

Remark  2.10  Corresponding  to  the  incident  field  <))o^\  for  |m|  >  n,  the  scat¬ 
tered  field  in  (21)  is  actually  dependent  on  n,  and  could  have  been  denoted 
there  by  In  other  words,  and  are  not  the  same  in  general, 

because  they  are  the  scattered  fields  corresponding  to  different  forward  models:  qn 
and  qn-i-  On  the  other  hand,  since  both  the  forward  models  look  the  same  as  the 
original  scatterer  q,  the  two  scattered  fields  and  have  the  same 

value  on  the  boundary  of  D{kw),  which  is  the  prescribed  scattering  data  (18). 
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2.5  A  Continuous  Version  of  the  Propagation  Number 

The  propagation  number  n  of  the  total  field 

(28) 

measures  the  depth  of  the  illuminated  area  A{n,  k):  the  smaller  it  is,  the  deeper 
the  total  field  penetrates  the  scatterer.  Since  A{n,k)  and  qn  depend  on  the 
discrete  variable  n,  it  is  analytically  improper  to  carry  out  perturbational  analysis 
of  A{n,  k)  and  on  n.  In  this  section,  we  generalize  the  propagation  number  n 
to  a  continuous  variable  rj. 

Given  a  real  number  6  >  0  (the  actual  value  of  b  is  qualitatively  unessential 
in  the  following  discussion),  consider  a  smooth  function  z/  :  defined  by 

the  formula 

'0,  7]  <  -b, 

1/(77)  =  •<  ^  >  0,  (29) 

monotone,  rj  €  0]. 

For  77  >  0,  we  denote  by  E{t])  the  scattering  experiment  with  the  weighted 
incident  fields 

=  i/{\Tn\  -  Tj)  ■  (30) 

for  all  integer  m.  We  further  denote  by  the  observable  part  of  q  to  the  exper¬ 
iment  E{tj)  (see  Remark  2.8),  and  by  the  corresponding  scattered  field,  so 

that 

= -iW  Gi(x,  {)■?,«)■[«!,"■'>«) (31) 

JD(w) 

It  is  easy  to  observe  the  following  facts. 

1.  If  77  >  6,  the  incident  fields  actually  used  in  the  experiment  E{q)  are  not 
all  but  only  those  whose  propagation  numbers  m  are  greater,  in  absolute  value, 
than  77  —  6. 

2.  If  77  =  N  is  substantially  greater  than  fco?,  every  incident  field  used  in  the 
experiment  E{rj)  can  only  penetrate  a  thin  layer  of  the  scatterer.  This  means 
that  for  all  |m|  >  N  —  6  the  Lippmann- Schwinger  equation  (31)  can  be  replaced 
by  the  Born  approximation  (see  Section  2.2) 

= -p  /  _  (32) 

JD{w) 

The  observable  forward  model  can  be  obtained  by  solving  these  linear  equa¬ 
tions  for  given  scattering  data 

{  |m|  >  N  -  6  }.  (33) 
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3.  At  t;  =  0,  incident  fields  of  all  propagation  numbers  are  used  in  the 
experiment,  and  they  are  equally  weighted  by  the  factor  v{\m\)  =  1.  Therefore, 
the  forward  model  qo  is  identical  to  qo,  both  being  the  observable  part  of  q 
corresponding  to  the  same  scattering  experiment  with  all  the  incident  fields 

{  for  all  integer  m  }.  (34) 

Remark  2.11  Clearly,  for  a  small  number  S,  the  scattering  experiment  E{t}  —  8) 
is  a  small  perturbation  to  the  experiment  Eirj').  More  specifically,  it  follows  from 
(SO)  that  the  incident  field  depends  on  rj  smoothly,  and  so  do  the  scattered 
and  total  fields,  due  to  the  well-posedness  of  the  forward  scattering  problem.  In 
the  reminder  of  the  paper,  we  assume  that  qr,  also  depends  on  q  smoothly. 


3  Layer-stripping  via  Skin  Effect 

The  aim  of  inverse  scattering  is  to  determine  9„|„=o  (see  Section  2.4,  and  Remark 
2.7  therein),  or  equivalently,  to  determine  qr,\r,=o  (see  Section  2.5).  The  two  for¬ 
ward  models  q,,  and  qn  depend  on  the  two  different  parameters  n  and  q,  one  being 
integer  and  therefore  discrete,  the  other  real  number  and  therefore  continuous. 
Accordingly,  the  iterative  inversion  method  to  be  presented  in  this  section  has 
two  versions:  one  determines  qo  via  perturbational  analysis  on  the  continuous 
parameter  q,  the  other  determines  on  the  discrete  parameter  n.  Their  de¬ 
scriptions  are  almost  identical.  Their  numerical  performances  are  similar.  Both 
versions  can  be  interpreted  as  layer-stripping  the  scatterer. 

3.1  Recursive  Linearization  on  tj 

The  nonlinear  inverse  problem  can  be  recursively  linearized  via  a  standard  per¬ 
turbational  analysis,  also  known  as  the  method  of  continuation,  on  q.  It  finds 
the  forward  model  qr,  first  for  a  large  q  where  the  nonlinear  equations  governing 
q,,  are  essentially  linear.  The  method  then  determines  q,,-s  for  a  small  number  6. 
The  nonlinear  equations  governing  qr,-s  can  be  linearized  by  those  governing  qr,. 
This  process  is  repeated,  and  finally,  is  obtained.  We  describe  the  inversion 
method  in  the  following  four  steps. 

Step  1  (Initialization).  Choose  N  sufficiently  greater  than  kw,  arid  divide 


vi  V2  m 

Figure  3:  Discretization  of  q. 
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the  interval  [0,n]  into  M  sufficiently  small  chunks  with  the  grid 


{  },  (35) 

where  t/o  =  0,  »/m  =  N,  and  rji  <  rji^i  for  0  <  i  <  M.  We  intend  to  obtain  qr,  at 
V  ~  VMt  ■  1  Vo- 


Step  2  (Born  Approximation).  At  7/  =  t/m  =  N,  obtain  the  observable  forward 
model  by  solving  the  linear  equations  (32). 

Remark  3.1  The  purpose  of  having  a  sufficiently  large  N  is  to  attain  the  skin 
effect  and  therefore  sufficiently  weak  scattering  (see  Section  2.2  for  details)  in 
order  to  make  the  Born  approximation  (32)  more  accurate.  It  turns  out  that  the 
Bom  approximation  needs  not  be  very  accurate  in  order  to  successfully  implement 
Step  2.  In  numerical  computation,  N  usually  can  he  chosen  only  slightly  greater 
than  kw. 

Remark  3.2  In  practice,  the  measurements  of  the  scattered  fields  (33)  on  the 
boundary  of  D{w)  are  not  available  for  |ml  ^  kw,  for  they  correspond  to  evanes¬ 
cent  waves;  therefore,  the  linear  equations  (32)  are  solved  for  finite  number  ofm. 
In  numerical  implementations,  only  two  linear  equations  with  m  =  ±N  are  solved 
to  produce  an  approximate  forward  model 


Step  3  (Linearization).  Now,  suppose  that  qij  has  already  been  obtained 
where  t)  >  0  is  a  point  on  the  grid  (35).  We  wish  to  determine  q^j  where  rj  is  the 

0  q  fj 

I - 1 - 1 - 1 - 1 - ! - 1 - i - 1 - 1 - 1 - 1 - 1 - \ - 1  I - 1 - ► 

Figure  4:  Update  from  rj  to  q. 

grid-point  immediately  to  the  left  of  q.  By  definition,  is  the  forward  model 
observable  to  the  scattering  experiment  E{q)  with  the  incident  fields  (30);  there¬ 
fore,  the  corresponding  scattered  field  is  connected  to  9,,  by  the  Lippmann- 
Schwinger  equation  (see  (31)) 

,  Gi(l,  ■  (36) 

JD(-a) 

In  (36),  %  and  are  unknown  inside  D{w),  whereas  u^’"’’’^(x)|r6aD(cj)  is  given 
as  the  scattering  data. 

3.1.  On  the  other  hand,  in  the  same  experiment  E{q),  the  forward  produces  a 
scattered  field,  denoted  here  by  corresponding  to  the  incident  field 
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The  scattered  field  satisfies  the  Lippmann- Schwinger  equation 

5(--’’)(x)  =  -k^  /  Gk{x, (37) 

JD{w) 

Since  qf,  is  known,  the  forward  problem  (37)  is  solved  for  the  scattered  field; 
therefore,  is  known  inside  as  well  as  on  the  boundary  of  D{w)^  for  every 

m  such  that  \m\  >  q  —  b. 

3.2.  By  assumption,  q  —  q  is  small  and  q^,  depends  on  q  smoothly  (see  Remark 
2.11);  therefore  the  perturbations 

^9  ~  Qn  ~  9i?i  (38) 

are  also  small.  For  every  m  such  that  \m\  >  q  -  b,  subtracting  (36)  from  (37), 
and  omitting  terms  quadratic  in  6q  and  we  obtain  an  equation  linking  8q 

linearly  to  8v^'^^: 

8v^^\x)  =  -k^  I  +  8q-[ut'^'>  +  (40) 

3.3.  Denoting  by  K,  Sm.  ■  c(D{w))  t-*  c{D{w))  the  two  compact  linear  operators 
defined  by  the  formulae 

K{w){x)  =  -k^  f  Gk{x,^yqi^{^yw{^)-d^,  (41) 

S„{w){x)  =  +  (42) 

for  arbitrary  w  €  c{D{w)),  we  rewrite  (40)  as 

8v^^^  =  K{8v^”^^)  +  Sm{8qy  (43) 

Since  the  forward  scattering  problem  is  well-posed,  the  operator  I  -  K  is  always 
invertible;  thus 

8v^^^  =  (I  -  K)-^  .  SmiSq).  (44) 

3.4.  Denote  by  P  :  c(D{w))  c(dD(w))  the  projection  operator,  and  by  ' 

c(D(w))  1-^  c(dD(w))  the  linear  operator  defined  by  the  formula 

Lm{w)  =  P  ■  {I  -  K)-^  .  S„,{w),  (45) 

for  arbitrary  w  €  c{D{w)).  Then  the  application  of  P  on  (44)  yields 

8v^'^\x)  =  Lra{8q),  (46) 
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for  all  X  €  dD{TX!),  and  all  m  such  that  |m|  >  Tj  —  b.  Since  the  linear  operator  can 
be  evaluated  and  therefore  is  known,  and  since  is  known  on  the  boundary 
dD{w)  (see  (36),  (37),  and  (39)),  the  linear  system  (46)  can  be  solved  as  a 
least-squares  problem  for  Sq.  Finally,  we  obtain  qr,  from  qij  and  6q  via  (38). 

Remark  3.3  In  reality,  is  not  available  for  |m|  kta  (see  Remark  3.2); 

therefore,  only  a  finite  numbers  of  the  linear  equations  (46)  need  to  be  solved.  In 
numerical  implementations,  the  linear  equations  with  77  —  6  <  |m|  <  N  are  solved 
to  produce  an  approximate  forward  model  q^ . 


Step  4  (Recursion).  The  process  of  linearization  described  in  Step  3  can 
obviously  be  iterated  in  two  directions.  First,  the  solution  %  is  accurate  to  the 
second  order  of  6q;  the  approximation  can  be  improved  by  replacing  q^  with  qr, 
and  repeating  the  computations  in  Step  3.  Secondly,  when  a  satisfactory  solution 
qr,  is  obtained  at  7/  =  77^,  it  is  used  in  Step  3  to  determine  the  next  forward  model 
qr,  aX  r)  =  rji-i  until  qo  is  found. 

Remark  3.4  It  turns  out  that  only  a  low  accuracy  of  the  solution  q,,  is  required 
in  numerical  computation  to  maintain  convergence  of  the  recursion  on  q.  In 
practice,  the  iteration  discussed  in  Step  4  to  improve  q^  is  often  unnecessary. 

3.2  Recursive  Linearization  on  n 

In  this  subsection,  we  describe  another  recursive  process,  almost  identical  to  that 
presented  in  the  preceding  subsection,  which  successively  determines  the  forward 
models 

qn,  71  =  N,N  -  1,...,1,0.  (47) 

This  linearization  process  is  accomplished  on  the  discrete  propagation  number 
n.  Its  procedures  are  parallel  to  those  of  the  preceding  subsection.  We  briefly 
present  them  here  for  completeness  and  for  the  reason  that  the  numerical  results 
given  in  Section  4  axe  produced  using  this  version  of  the  linearization  method. 

Step  1  (Initialization).  Choose  N  sufficiently  greater  than  kw,  so  that  the 
skin  effect  is  attained  (in  practice,  N  is  usually  chosen  just  slightly  greater  than 
kw,  see  Remark  3.1). 

Step  2  (Born  Approximation).  Obtain  the  observable  forward  model  q^  by 
solving  the  linear  equation  (10)  with  m  =  ±N. 

Step  3  (Linearization).  Now,  suppose  that  qn+i  has  already  been  obtained 
where  ti  <  N  is  a  positive  integer.  We  wish  to  determine  qn  which  is  the  for¬ 
ward  model  observable  to  the  scattering  experiment  E{n),  and  which  satisfies 
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the  equations  (see  (21)) 

-e  I  (48) 

JA(n,k) 

for  all  m  such  that  |m|  >  n.  In  (48),  is  unknown  inside  D{w),  but  is  given 
on  the  boundary  of  D{w)  as  the  scattering  data  (see  (23)). 

3.1.  In  the  same  scattering  experiment  E{n),  the  forward  model  ^n+i,  which  is 
known,  produces  the  scattered  field,  denoted  here  by  corresponding  to  the 
incident  field  Therefore,  satisfies  the  equation 

-P  /  (49) 

and  is  obtained  by  solving  the  forward  problem  (49). 

3.2.  Assuming  that  is  close  to  (see  Remark  3.6),  we  observe  that  the 
perturbations 

Sq  =  qn+i  -  qn,  (50) 

(51) 

are  small.  For  every  |m|  >  n,  subtracting  (48)  from  (49),  and  omitting  terms 
quadratic  in  6q  and  we  obtain  an  equation  linking  8q  linearly  to 

„  Gt(i,f){?n+i-«i/’<">  +  (52) 

JA{n,k) 

3.3.  Denoting  by  K,  Sm  '  c(A(n,A:))  i— >  c{A{n,k))  the  two  compact  linear  oper¬ 


ators  defined  by  the  formulae 

K{w){x)  =  -P/  Gfc(x,0-9n+i(0-M0-<^^>  (53) 

JA(n,k) 

S„(w)(x)  =  -el  ^  Gt(i,0-[4”>«)  +  ^<"l(0)-to(0<,  (54) 

J  A(n,k) 

for  arbitrary  w  G  c(A{n,  k)),  we  rewrite  (52)  as 

8tP^^^  =  K{8i,^^^)  +  S„,i8q).  (55) 

Since  the  forward  scattering  problem  is  well-posed,  the  operator  I  —  K  is  always 
invertible;  thus 

8x1;^^^  =  {I  -  K)-^  ■  Sm{8q).  (56) 

3.4.  Denote  by  P  :  c{A{n,k))  i— »  c{dD{zu))  the  projection  operator,  and  by 
Lm  '■  c{A{n,  k))  c{dD{w))  the  linear  operator  defined  by  the  formula 

L^{w){x)  =  P-{I-  K)-^  •  5„.(u;),  (57) 
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(58) 


for  arbitrary  w  e  c(A(n,  k)).  Then  the  application  of  P  on  (56)  yields 

=  L„(6q), 

for  all  X  €  dD{w),  and  all  m  such  that  |m|  >  n.  Since  the  linear  operator  can 
be  evaluated  and  therefore  is  known,  and  since  is  known  on  the  boundary 

dD{w)  (see  (48),  (49),  and  (12)),  the  linear  system  (58)  can  be  solved  as  a 
least-squares  problem  for  6q.  Finally,  we  obtain  qn  from  g„+i  and  6q  via  (51). 

Step  4  (Recursion).  Repeat  the  linearization  procedures  in  Step  3  till  n  =  0 
where  qo  is  obtained. 

Remark  3.5  In  reality,  is  not  available  for  |m|  ^  kw  (see  Remark  3.2); 

therefore,  only  finite  numbers  of  the  linear  equations  (58)  need  to  be  solved.  In 
numerical  implementation,  the  equations  with  n  <  |m|  <  N  are  solved  to  produce 
an  approximate  forward  model  qn. 

Remark  3.6  The  perturbation  Sq  of  (50),  unlike  that  of  (38),  can  not  be  made 
arbitrarily  small,  since  n  in  (50)  is  discrete.  In  fact,  our  numerical  experiment 
shows  that  Sq  may  not  be  small  at  all.  This,  in  principle,  makes  the  perturbational 
analysis  invalid  and  may  lead  to  divergence  of  the  recursive  linearization.  It  turns 
out  that  Sq  is  small  in  A{n-\~\,  k),  the  illuminated  area  by  the  experiment  E{n+1); 
therefore,  the  perturbational  analysis  can  be  carried  out  there  to  linearize  the 
relationship  between  ^9U(n+i,fc)  and  However,  Sq  usually  is  not  small  in 

T{n,k)  =  A{n,k)\A{n  +  l,k),  (59) 

which  is  the  area  illuminated  by  E{n)  but  not  by  .S(n-t-l),  and  which  we  refer  to  as 
the  transition  zone.  Fortunately,  this  layer  is  dimly  illuminated  by  the  experiment 
E{n),  and  across  it  the  skin  effect  occurs.  Therefore,  in  the  transition  zone,  the 
Bom  approximation  is  valid  (see  Section  2.2)  and  it  linearizes  the  relationship 
between  ^|T(n,Jt)  and  the  scattered  fields  of  (48);  consequently,  the  relationship 
between  Sq\x{n,k)  =  9|r(n,A:)  and  is  also  essentially  linear.  In  summary,  in 

spite  of  the  fact  that  6q  may  not  be  small  in  A{n,  k),  it  is  essentially  linearly 
connected  to  as  is  formulated  in  (52)  and  finally  in  (58). 

3.3  Recursive  Linearization  with  an  Initial  Guess 

Neither  of  the  inversion  procedures  described  in  Sections  3.1  and  3.2  requires  an 
initial  guess  of  the  scatterer  to  start.  But  an  initial  guess  can  be  incorporated  into 
the  inversion  procedures  and  proves  beneficial  to  the  reconstruction  (see  Example 
3,  Section  4)  if  it  is  not  a  bad  approximation  to  the  scatterer  to  be  recovered. 
Usually,  an  initial  guess  is  produced  with  a  method  of  inversion  at  a  relatively 
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low  frequency  k,  where  low-frequency  components  of  the  scatterer  are  recovered 
(see  [7]  for  more  details  and  numerical  examples).  The  inversion  method  used  to 
obtain  an  initial  guess  could  be  the  one  presented  in  [7]  or  described  here  in  this 
paper. 

The  introduction  of  the  initicil  guess  to  the  two  inversion  procedures  of  Sec¬ 
tions  3.1  and  3.2  is  straightforward.  We  take  the  latter  as  an  example  to  describe 
it.  The  only  modification  is  made  in  Step  2  of  Section  3.2  where  is  written  as 

9n  =  ^initial  +  ^<1-  (60) 

Consequently,  6q  can  be  obtained  by  solving  the  linear  equation  (10)  with  m  = 
±N.  With  5n  known,  we  proceed  to  Step  3  without  any  further  alterations  to  the 
algorithm. 

4  Numerical  Experiments  with  Radially  Sym¬ 
metric  Scatterers 

Both  versions  of  the  recursive  linearization  method  (see  Sections  3.1  and  3.2) 
were  implemented  numerically.  But  since  their  performances  are  very  similar, 
we  only  present  numerical  results  of  the  second  version  based  on  the  discrete 
propagation  number  n. 

Although  the  inverse  scattering  problem  and  its  linearization  are  described 
in  this  paper  in  terms  of  the  Lippmann-Schwinger  integral  equations,  models 
based  on  other  equations  exist,  for  example,  the  Riccati  equations  (see  [6],  [7]). 
Since  an  accurate  and  reliable  code  solving  the  integral  equations  is  not  presently 
available  to  the  author,  the  actual  equations  employed  are  the  Riccati  equations 
whose  Born  approximation  (required  in  Step  2  of  Section  3.2)  and  perturbational 
analysis  (all  the  procedures  in  Step  3  of  Section  3.2)  are  discussed  in  great  details 
and  implemented  numerically  in  [7]. 

For  reconstruction  of  a  general  two-dimensional  scatterer  of  N  x  N  wave¬ 
lengths,  it  turns  out  that  the  algorithm  based  on  the  integral  equations  requires 
0{N^)  operations,  whereas  the  algorithm  based  on  the  Riccati  equations  requires 
0{N^)  operations,  an  approach  too  expensive  to  be  used  in  solving  a  problem 
with  N  =  10.  In  this  section,  therefore,  we  only  present  numerical  results  of 
inversion  of  cylindrically  symmetric  scatterers,  for  which  the  algorithm  requires 
0{N^)  operations.  The  numerical  experiments  for  general  scatterers  will  be  pre¬ 
sented  later  when  the  code  for  the  integral  equations  is  acquired 

4.1  The  Forward  Problem 

As  is  mentioned  in  the  beginning  of  Section  4,  the  Riccati  equations,  instead  of  the 
Lippmann-Schwinger  equations,  are  used  in  the  forward  modeling.  A  FORTRAN 
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code  was  written  to  solve  the  forward  problem  required  in  Step  3.1  of  Section 
3.2;  see  [7]  for  details  of  the  numerical  treatment  of  the  Riccati  equations.  For  a 
prescribed  scatterer  q,  the  forward  problem  is  also  solved  to  produce  the  synthetic 
scattering  data  (see  (23)).  The  numerical  solutions  in  both  cases  are  accurate 
at  least  to  four  digits.  For  a  cylindrically  symmetric  scatterer,  the  forward  solve 
requires  0(N‘^)  operations. 


4.2  The  Least-squsires  Problem 

Since  the  ill-posedness  of  the  original  nonlinear  inverse  problem  is  inherited  by  the 
linearized  problem  (58),  it  must  be  solved  as  a  least-squares  problem.  This  is  easy 
because  it  is  a  linear  problem,  for  which  the  standard  conjugate  gradient  method 
was  used  (see  [7]  for  more  details).  As  with  solving  other  least-squares  problems, 
the  conjugate  gradient  iterations  must  be  carefully  measured  and  terminated  for 
the  stability  of  the  solution.  Here,  the  stop  condition  is  on  the  magnitude  of 
the  residual:  the  iterative  procedure  is  terminated  and  the  solution  is  obtained 
when  the  norm  of  the  residual  is  comparable  to,  or  not  less  than,  the  pre¬ 
cision  in  which  the  forward  problem  is  solved  or  the  scattering  measurement  is 
collected.  This  is  the  only  place  in  the  inversion  method  where  “regularization 
of  ill-posedness”  is  numerically  imposed. 

4.3  Numerical  Results 

Several  technical  details  regarding  the  presentation  of  the  numerical  results  re¬ 
quire  explanations. 

1.  The  radius  zu  of  the  disk  scatterer  is  chosen  as  w  =  n  for  convenience, 
so  that  the  number  of  wavelengths  across  the  diameter  of  the  disk  is  the  wave 
number  k. 

2.  In  plots  showing  the  results  of  reconstructions,  the  prescribed  scatterer 
functions  axe  depicted  in  solid  curves,  whereas  the  numerically  reconstructed 
scatterers  are  represented  by  dotted  curves. 

3.  The  error  of  a  reconstruction  is  defined  by  the  formula 

^2  =  ^  [9(^)  -  qoY-r-dr^ '  ,  (61) 

where  is  the  reconstructed  scatterer  obtained  at  the  end  of  the  recursive  algo¬ 
rithm  (see  Section  3.2). 

4.  Scattering  data  of  various  accuracies  are  used  in  the  inversion  to  examine 
the  stability  of  the  algorithm.  Scattering  data  of  low  precision  is  obtained  from 
data  of  high  precision  by  truncation.  We  say  that  the  data  is  accurate  to  D  digits 
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if  each  number  in  the  high  precision  data  is  truncated  to  D  digits.  For  instance, 
truncating  the  number  0.12345  to  three  digits  yields  0.123. 

Example  1.  Reconstruct  a  scatterer  defined  by  the  formula 

q{r)  =  0.1  •  [1  -  cos(14r)],  (62) 

inside  D{'k).  Four  inversions  were  made  aX  k  =  8.6  using  scattering  data  accurate 
to  Z)  =  1,  2,  3,  4  digits;  see  Table  1  for  the  error  of  the  reconstructions.  Figure 
5  shows  the  reconstructed  scatterer  against  the  exact,  for  Z>  =  1,  2,  3,  4.  With 
Z)  =  3,  the  process  of  layer-stripping  is  evidently  shown  in  Figures  6 — 8  where 
the  reconstructions  of  qn  of  (47)  are  plotted. 


D 

1 

2 

3 

4 

62 

0.7450E-1 

0.1719E-1 

0.648E-2 

0.598E-2 

Table  1:  Errors  of  Reconstructions  with  D  =  1,2, 3,4,  Example  1. 


Example  2.  Reconstruct  a  scatterer  defined  by  the  formula 

9(r)  =  0.26  •  [0.55cos(4r)  -  0.44sin(8r)  -1-  0.25sin(12r)  +  0.14cos(16r)],  (63) 

inside  Z?(7r).  Four  inversions  were  made  at  frequency  k  =  9.3  with  data  accu¬ 
rate  to,  respectively,  I>  =  1,  2,  3,  4  digits;  see  Table  2  for  the  L'^  error  of  the 
reconstructions.  Four  more  inversions  were  also  made  separately  at  frequencies 


D 

1 

2 

3 

4 

62 

0.181 

0.421E-1 

0.309E-1 

0.198E-1 

Table  2:  Errors  of  Reconstructions  sX  k  =  9.3,  D  =  1,2, 3,4,  Example  2. 

k  =  2.0,  4.1,  6.2,  9.3  with  scattering  data  accurate  to  four  digits.  Table  3  shows 
the  error  of  the  reconstructions;  Figures  9  and  10  show  plots  of  the  recon¬ 
structions  and  the  error  distributions.  According  to  Heisenberg’s  uncertainty 


k 

2.0 

4.1 

6.2 

9.3 

62 

0.691 

0.332 

0.155 

0.172E-1 

Table  3:  Errors  of  Reconstructions  at  Various  Values  of  k,  Example  2. 
principle  (see  [7]),  only  those  Fourier  modes  of  the  scatterer  q  whose  frequencies 
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are  lower  than  or  about  2k  can  be  recovered  in  a  reconstruction  at  wave  number 
k.  Therefore,  at  fc  =  2.0,  only  the  first  term  of  (63)  is  recovered;  at  k  =  4.1, 
roughly  the  first  two  terms  of  (63)  are  recovered;  zX  k  =  6.2,  roughly  the  first 
three  terms  of  (63)  are  recovered;  and  ai  k  =  9.3,  all  four  terms  are  essentially 
recovered.  In  fact,  if  we  regard  the  first  term  alone,  the  first  two  terms  alone, 
and  the  first  three  terms  alone  as  approximations  to  the  scatterer  function  (63), 
then  they  approximate  (63)  with  errors  0.701,  0.339,  and  0.170,  which  closely 
match  the  errors  of  inversions  given  in  Table  3. 

At  A:  =  4.1,  6.2,  the  reconstructions  of  the  intermediate  forward  models  qn  of 
(47)  in  the  recursive  linearization  process  are  shown  in  Figures  11  and  12.  The 
corresponding  plots  at  A:  =  2.0,  9.3  are  not  shown  here  because  they  are  similar 
to  those  shown  in  Figures  11  and  12. 


Example  3.  Reconstruct  a  scatterer  defined  by  the  formula 


q{r)  =  0.2  •  [sin(4r)  +  cy  sin(8r)  +  C2-sin(14r)],  (64) 


with  Cl  =  —0.846153846,  C2  =  0.153846154.  The  scatterer  was  reconstructed  at 
k  =  9.7.  Figure  13  shows  plots  of  the  reconstruction  and  the  error  distribution. 
Four  more  inversions  were  made  with  four  different  initial  guesses  qiniti,  qinit2, 
qinita,  and  qinit4  defined  by  the  formulae 


qinitiir) 

?tmt2(^) 

9tnit3(7’) 

qinitiij’') 


0.2  •  [sin(4r)  —  ci-  sin(8r)],  (65) 

0.2  •  [sin(4r)  —  ci-  sin(8r)  +  C2-  sin(14r)  +  0.1-  sin(13r)],  (66) 

0.2  •  [sin(4r)  —  Ci-  sin(8r)  +  C2-  sin(14r)  +  O.T  sin(18r)].  (67) 

0.2  •  [sin(4r)  —  Ci-  sin(8r)  +  C2-  sin(14r)  +  0.1-  sin(22r)].  (68) 


Table  4  shows  the  errors  of  the  initial  guesses  and  of  the  final  reconstructions. 
Again  according  to  Heisenberg’s  uncertainty  principle  (see  [7]),  the  highest  fre¬ 
quency  of  the  Fourier  modes  of  the  scatterer  that  can  be  recovered  is  about 
2A;  =  2  X  9.7  =  19.4;  therefore,  the  mode  sin(13r)  or  sin(18r)  in  the  initial  guess 
is  recognizable  in  the  scattering  experiment  and  can  be  eliminated  from  the  initial 
guess  in  the  reconstruction,  whereas  the  mode  sin(22r)  in  qinit4  is  not  observable 
to  the  scattering  experiment  and  can  not  be  eliminated  from  the  initial  guess  in 
the  reconstruction;  see  Table  4. 


Init.  Guess 

None 

Qiniti 

9init2 

^initS 

62  (Init.) 

1.0 

0.114 

0.744E-1 

62  (Final) 

0.133E-1 

0.529E-2 

0.262E-2 

0.138E-1 

0.673E-1 

Table  4:  Errors  of  Reconstructions  with  Initial  Guesses,  Example  3. 
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Example  4.  Reconstruct  a  scatterer  defined  by  the  formula 


_  J  0-2  •  [sin(4r)  +  cy  sin(8r)  +  C2-  sin(14r)],  r  <  2.2707,  .  . 

^  \  0  2.2707  <  r  <  TT, 

with  Cl  =  —0.846153846,  C2  =  0.153846154.  The  scatterer  function  is  identical 
to  that  of  Example  3  for  r  <  2.2707,  and  there  is  a  discontinuity  across  the  circle 
r  =  2.2707.  The  inversion  was  made  at  A:  =  9.7.  As  is  expected,  the  Gibbs’ 
phenomenon  appears  in  the  reconstructed  scatterer  near  the  discontinuity;  see 
Figure  14  for  plots  of  the  reconstruction  and  the  error  distribution.  The  recovery 
of  the  intermediate  forward  models  of  (47)  are  shown  in  Figures  15 — 17. 


5  Discussions 

We  have  presented  a  new  inversion  method  and  its  preliminary  numerical  results 
in  the  special  case  of  cylindrically  symmetric  scatterers.  In  this  section,  we  make 
several  technical  remarks  on  the  numerical  experiments  provided  in  Section  4.3, 
and  discuss  possible  directions  of  future  investigation  of  the  method. 

1.  Numerical  implementations  for  general  scatterers  in  two  dimensions  will 
be  presented  on  a  later  date  when  a  code  solving  the  two-dimensional  Lippmann- 
Schwinger  equation  is  available.  The  extension  of  the  inversion  method  to  the 
problem  of  electrical  impedance  imaging,  though  straightforward,  constitutes  a 
sepaxate  topic,  and  will  be  investigated  and  presented  separately.  In  this  case,  a 
disk  containing  the  inhomogeneity  is  probed  with  electrical  potential  of  the  form 

<l>o{r,e)  =  r”^  m  >  0.  (70) 

Therefore,  every  probing  field  with  m  >  0  is  a  decaying  mode  in  any  interval 
[0,  A].  This  means  that  the  number  of  parameters  which  can  be  recovered,  inde¬ 
pendent  of  the  size  of  the  scatterer,  is  limited.  In  many  applications,  even  such 
a  limited  resolution  is  difficult  to  acquire,  and,  if  attainable,  would  be  useful. 

2.  The  skin  effect  and  the  illuminated  areas  of  various  depth  can  be  achieved 
in  three  dimensions  and  in  geometries  other  than  that  of  the  disk  used  here  in 
this  paper.  For  example,  in  underwater  environment  where  the  transducers  are 
restricted  on  the  surface  of  water,  we  may  use  the  incident  waves  of  the  form 

^o(x,y,^)  =  e‘■”‘"+•■""•e-P^  (71) 

to  attain  the  skin  effect,  where  is  the  depth  of  the  ocean,  m,  n  are  real  numbers 
such  that  +  ti^  >  so  that  p  =  y/m^  +  —  k'^  is  a  positive  number.  There¬ 

fore,  p  >  0  is  a  parameter  which  controls  the  depth  of  penetration  of  the  incident 
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field.  Recursive  linearization  can  be  carried  out  on  the  continuous  parameter  p 
for  decreasing  p  till  p  =  0.  Then  the  recursion  may  be  on  the  real  parameter  I  of 
the  incident  field 

=  (72) 

for  increasing  I  in  the  interval  [0,  k]. 

3.  The  stability  of  the  method  is  demonstrated  by  reconstructions  with  noise- 
corrupted  scattering  data  in  Examples  1  and  2.  More  stability  tests  were  made, 
the  results  being  similar  to  those  presented  here.  The  inversion  method  does 
not  require  the  solution  of  unstable  initial  value  problems  which  the  standard 
layer-stripping  algorithm  must  solve  in  order  to  recover  the  scatterer  layer  by 
layer.  When  the  ill-posed  and  nonlinear  inverse  problem  is  recursively  linearized, 
the  ill-posedness  is  passed  to  each  linear  problem  which  we  know  very  well  how 
to  solve  to  maintain  stability.  Thus,  the  issue  of  stability  is  trivial  here  so  far  as 
we  do  not  insist  on  recovering  modes  of  the  scatterer  that  are  not  observable  to 
our  scattering  experiment.  Conversely,  if  somehow  during  the  recursive  precess 
the  intermediate  result  is  polluted  with  such  unobservable  modes,  usually  high- 
frequency  modes,  there  is  no  mechanism  of  the  algorithm  that  can  remove  them, 
as  is  shown  in  Example  3. 

4.  The  algorithm  starts  by  illuminating  the  scatterer  with  an  incident  wave 
of  high  propagation  number  that  penetrates  only  a  thin  layer  of  the  scatterer. 
There,  the  measurement  of  the  scattered  field  is  weak,  and  therefore  is  prone  to 
contaminations  of  noise.  In  principle,  the  propagation  number  m  is  required  to  be 
considerably  greater  than  kw  (see  Section  2.2)  in  order  to  make  the  relationship 
between  the  measurement  and  the  scatterer  essentially  linear.  In  practice,  fortu¬ 
nately,  only  a  roughly  linear  relationship  is  needed  to  produce  an  approximate  ob¬ 
servable  forward  model  in  the  first  layer.  More  precisely,  the  highest  propagation 
number  m  can  be  chosen  only  slightly  greater  than  kw  such  that  the  measure¬ 
ment  of  the  far-field  can  be  10% — 30%  as  strong  as  the  strongest  measurement 
of  all  propagation  numbers.  Thus,  the  so-called  weak  scattering  required  theo¬ 
retically  to  attain  ideal  linearization  does  not  practically  cause  problem  to  the 
noise-to-signal  ratio:  if  the  strongest  scattering  can  be  measured  accurately,  so 
can  the  weakest. 
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Figure  5:  Reconstructions  &t  k  =  8.6,  D  =  1,2, 3, 4,  Example  1. 
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Figure  6:  Reconstructions  of  qn  aX  k  =  8.6,  n  =  29, 28, . . . ,  20,  Example  1. 
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Figure  9:  Reconstructions  at  A:  =  2.0,  4.1,  6.2,  9.3,  Example  2. 


Figure  10:  Errors  of  Reconstructions  at  A;  =  2.0,  4.1,  6.2,  9.3,  Example  2. 
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Figure  11:  Reconstructions  of  at  A:  =  4.1,  n  =  15, 14, . . . ,  0,  Example  2. 
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Figure  16:  Reconstructions  of  qr,  &t  k  =  9.7,  n  =  19, 18, . . . ,  10,  Example  4 
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